FORMING PLANETESIMALS BY GRAVITATIONAL INSTABILITY. II. HOW DUST SETTLES TO ITS MARGINALLY STABLE STATE

Dust can settle quickly in gaseous protoplanetary disks. In a passive (non-turbulent) nebula, a particle's vertical height z above the midplane obeys

$$\begin{equation} \ddot{z} = -\dot{z}/t_{\rm stop}- \Omega _{\rm K}^2 z, \end{equation} \tag{ 1 }$$

where the first term on the right-hand side accounts for gas drag and the second term accounts for stellar gravity when z ≪ r, the cylindrical radius. Here, Ω_K is the Keplerian orbital frequency and

$$\begin{equation} t_{\rm stop}\equiv \frac{m v_{\rm rel}}{F_{\rm D}} \end{equation} \tag{ 2 }$$

is the momentum stopping time of a particle of mass m moving at speed v_rel relative to gas. Expressions for the drag force F_D can be found in Adachi et al. (1976) and Weidenschilling (1977). We are interested in small, well-coupled particles having non-zero stopping times much shorter than the dynamical time: 0 < τ_s ≡ Ω_Kt_stop ≪ 1. Spherical particles of radius s and internal density ρ_s that experience Epstein drag (F_D ∝ s²v_rel so that t_stop does not depend on v_rel) settle to the midplane at terminal velocity −Ω²_Kzt_stop in a time

$$\begin{equation} t_{\rm settle} \sim \frac{1}{\Omega _{\rm K}\tau _{\rm s}} \sim 10^3 \left(\frac{0.1 {\rm \,cm}}{s} \right) \left(\frac{1 {\rm \,g}{\rm \,cm}^{-3}}{\rho _{\rm s}} \right) \left(\frac{F}{1} \right) {\rm \,yr}\,. \end{equation} \tag{ 3 }$$

For this and all other numerical evaluations in this paper, we use a background disk that is F times more massive than the minimum-mass nebula of Chiang & Youdin (2010, hereafter CY10; see Appendix A). For such a disk, t_settle is nearly independent of stellocentric distance. The assumption that particles are spherical may not be too bad because fractal aggregates of grains are expected to become more compact as they collide with one another (Dominik & Tielens 1997; Dullemond & Dominik 2005; Ormel et al. 2007).

For millimeter-sized particles, the settling time t_settle is much shorter than the disk lifetime, measured in Myr (Hillenbrand 2005; Hernández et al. 2008). By comparison, micron-sized and smaller particles stay suspended at least one scale height above the midplane as long as the gas disk is present. To the extent that collections of particles of different sizes tend to place their mass at the upper end of the size distribution and their surface area at the lower end, we can expect most of the solid mass in disks to sediment out into a thin sublayer, leaving behind the smallest of grains to absorb incident starlight in a flared disk atmosphere. On the whole this picture is consistent with observed spectral energy distributions of T Tauri disks, although some models hint that large grains might remain lofted up in a disk two gas scale heights thick (D'Alessio et al. 2006). Settling can only proceed when and where disk turbulence dies, in regions where gas is insufficiently dense to sustain gravitoturbulence (Gammie 2001) and too poorly ionized to be magnetorotationally unstable (Gammie 1996). Our current understanding of disk turbulence easily admits such passive disk regions, located at disk midplanes where ionizing radiation from the star cannot penetrate (Perez-Becker & Chiang (2010). Turner et al. (2010) found in numerical simulations that even when disk surface layers were magnetorotationally unstable, grains at the midplane settled much as they would in a laminar flow.

In passive disk regions, how far down do dust particles settle? When particles are small enough not to be affected by aerodynamic streaming instabilities (SIs; e.g., Bai & Stone 2010a), we expect them to settle until the dust density gradient ∂ρ_d/∂z becomes so large and the consequent vertical shear in orbital velocity ∂v_ϕ/∂z becomes so strong that the sublayer is on the verge of overturning by a Kelvin–Helmholtz-type instability (KHI; Weidenschilling 1980). An order-of-magnitude estimate of the minimum layer thickness can be derived using the Richardson number

$$\begin{equation} Ri \equiv \frac{-(g/\rho)(d\rho /dz)}{(dv_\phi /dz)^2} \,, \end{equation} \tag{ 4 }$$

which if less than some critical value Ri_crit may signal that the layer is KH unstable (e.g., Drazin & Reid 2004). Here, g is the vertical gravitational acceleration and ρ = ρ_d + ρ_g is the total density of dust plus gas. In Lee et al. (2010, hereafter Paper I), we found that Ri_crit increases with Σ_d/Σ_g, the ratio of dust-to-gas surface densities, a.k.a. the bulk metallicity. For disks of bulk solar metallicity, we determined empirically that Ri_crit ≈ 0.2.

To translate the Richardson number (Equation (4)) into a critical dust layer thickness, first recognize that the orbital velocity v_ϕ depends on the local dust-to-gas ratio μ ≡ ρ_d/ρ_g according to

$$\begin{equation} v_\phi = \Omega _{\rm K}r \left(1 - \frac{\eta }{\mu +1} \right) \end{equation} \tag{ 5 }$$

in the inertial frame, where

$$\begin{eqnarray} \eta & \equiv & \frac{ -(1/\rho _{\rm g}) \partial P / \partial r }{2 \Omega _{\rm K}^2r} \approx \frac{1}{2} \left(\frac{c_{\rm s}}{\Omega _{\rm K}r} \right)^2 \approx \frac{1}{2} \left(\frac{H_{\rm g}}{r} \right)^2 \nonumber \\ & \approx & 8 \times 10^{-4} \left(\frac{r}{\rm AU} \right)^{4/7} \end{eqnarray} \tag{ 6 }$$

is a dimensionless measure of the strength of the background radial pressure gradient ∂P/∂r, with gas scale height H_g and sound speed c_s (e.g., Nakagawa et al. 1986). When ∂P/∂r < 0, pressure provides extra support against radial stellar gravity and so drives the gas to move on slower than Keplerian orbits. The orbital velocity depends on μ as in Equation (5) because dust-laden gas, weighed down by the extra inertia of solids, is accelerated less by the radial pressure gradient than is dust-free gas, and so must hew more closely to Keplerian rotation. Call the critical layer height Δz_Ri for which Ri = Ri_crit, and assume the midplane gas-to-dust ratio μ₀ ≳ 1 (above the layer μ ≪ 1). Then evaluating Equation (4) with the approximations g ≈ −Ω²_KΔz_Ri (no self-gravity), ρ⁻¹∂ρ/∂z ∼ −1/Δz_Ri, and ∂v_ϕ/∂z ∼ −ηΩ_Kr/Δz_Ri, we find

$$\begin{eqnarray} \Delta z_{\rm Ri} & \sim & Ri_{\rm crit}^{1/2} \eta r \sim \frac{1}{2} Ri_{\rm crit}^{1/2} \frac{H_{\rm g}}{r} H_{\rm g} \nonumber \\ & \sim & 5 \times 10^{-3} \left(\frac{Ri_{\rm crit}}{0.2} \right)^{1/2} \left(\frac{r}{\rm AU} \right)^{2/7} H_{\rm g} \,. \end{eqnarray} \tag{ 7 }$$

Equation (7) indicates that the dust layer could be quite thin, subtending on the order of 1% of the gas scale height. Is this thin enough for the dust to self-gravitate and hopefully fragment into planetesimals? One can compare the midplane density to the "Toomre density" required for the disk to undergo gravitational instability on the dynamical time Ω⁻¹_K (Safronov 1969; Goldreich & Ward 1973).⁵ As reviewed by CY10, the Toomre density is⁶

$$\begin{equation} \rho _{\rm Toomre} \approx \frac{M_{\ast }}{2\pi r^3} \approx 10^{-7} \left(\frac{r}{{\rm AU}} \right)^{-3} {\rm \,g}{\rm \,cm}^{-3} \end{equation} \tag{ 8 }$$

where the numerical evaluation is for a central stellar mass M_* equal to 1 M_☉. Now the actual midplane (subscript 0) density is

$$\begin{equation} \rho _0 = \rho _{\rm g0} + \rho _{\rm d0} = 2.7 \times 10^{-9} F (1 + \mu _0) \left(\frac{r}{\rm AU} \right)^{-39/14} {\rm \,g}{\rm \,cm}^{-3} \end{equation} \tag{ 9 }$$

which means the midplane dust-to-gas ratio must be

$$\begin{equation} \mu _{\rm 0, Toomre} \approx 34 \left(\frac{1}{F} \right) \left(\frac{M_{\ast }}{M_\odot } \right) \left(\frac{r}{\rm AU} \right)^{-3/14} \end{equation} \tag{ 10 }$$

for the midplane density to match the Toomre density. By comparison, in our crude model of a dust sublayer whose height above the midplane cannot be smaller than Δz_Ri, the midplane dust-to-gas ratio cannot exceed

$$\begin{equation} \mu _{\rm 0,Ri} \sim \frac{ \Sigma _{\rm d}/ (2\Delta z_{\rm Ri}) }{\rho _{\rm g0}} \sim 1 \left(\frac{\Sigma _{\rm d}/\Sigma _{\rm g}}{0.015} \right) \left(\frac{0.2}{Ri_{\rm crit}} \right)^{1/2} \left(\frac{r}{\rm AU} \right)^{-4/7} \end{equation} \tag{ 11 }$$

which is nominally smaller than μ_0,Toomre by more than an order of magnitude. Here the bulk (height-integrated) metallicity Σ_d/Σ_g is normalized to solar abundance (Lodders 2003), assuming all metals have condensed into grains.⁷

For many years, the fact that μ_0,Ri falls short of μ_0,Toomre was believed to rule out the formation of planetesimals by collective effects, self-gravitational or otherwise (e.g., Weidenschilling & Cuzzi 1993), but there are more ways to achieve the Toomre density than vertical settling. A dissipative form of gravitational instability can, in principle, collect particles radially into overdense rings even when self-gravity is weaker than stellar tidal forces (Ward 1976, 2000; Coradini et al. 1981; Youdin 2005; for a simple explanation of the instability, see the introduction of Goodman & Pindor 2000). It is not clear whether this instability, which operates over lengthscales and timescales longer than those characterizing the Toomre instability by at least a factor of (μ_0,Toomre/μ₀)², can compete with other effects that seek to rearrange dust and gas (e.g., Youdin 2005).

Another alternative is to invoke larger dust particles that are only marginally coupled to gas; these can clump by the aerodynamic SI (Youdin & Goodman 2005; Johansen et al. 2009; Bai & Stone 2010a, 2010c). In their three-dimensional (3D) numerical simulations, Johansen et al. (2009) reported that particles having τ_s = 0.1–0.4—corresponding to sizes of a few centimeters at r = 5 AU if F = 1, and larger sizes if F>1—concentrated so strongly by aerodynamic effects that planetesimals effectively hundreds of kilometers across coalesced within just a few orbits. To obtain this result, Johansen et al. (2009) initialized their simulations by placing the bulk of the disk's solid mass into particles approaching decimeters in size. Bai & Stone (2010a) greatly expanded the range of τ_s modeled and found similar results for their 3D simulations: in the highly turbulent states driven by the SI, instantaneous densities exceeded the Roche density⁸ when the disk's solids were all composed of particles having τ_s = 0.1–1 and the bulk height-integrated metallicity was about twice solar; see run R10Z3-3D in their Figure 5. For this same run, the time-averaged dust-to-gas ratio at the midplane was ∼12, a factor of a few less than the Toomre threshold; see their Figure 4 and compare with our Equation (10). By contrast, when half or more of the disk's solid mass had τ_s < 0.1, or when disks had smaller metallicities, their simulated densities fell short of the Roche and Toomre densities by more than an order of magnitude. Note that we are quoting from the 3D simulations of Bai & Stone (2010a).

Given how sensitive the SI is to the existence of marginally coupled particles (centimeter to meter sized for τ_s ∼ 0.1–1, r ∼ 1–30 AU, and order unity F), whether enough such particles actually exist in protoplanetary disks for the SI to play a dominant role in planetesimal formation remains an open and delicate issue. Appeal is often made to observed spectral energy distributions and images of T Tauri disks at centimeter wavelengths; these suggest that much of the solid mass is in millimeter-to-centimeter sized particles (e.g., D'Alessio et al. 2001; Wilner et al. 2005). Larger-sized particles are plausibly also present but are not inferred for want of data probing the disk at longer wavelengths. One problem concerns how quickly τ_s ≳ 0.1 particles can be grown, and how they can survive orbital decay by gas drag. In the 3D simulations of Bai & Stone (2010a), the SI clumped particles strongly enough for self-gravity to be significant when τ_s ≳ 0.1 particles comprised more than half of the disk's solid mass. It is unclear whether particle–particle sticking can build up such a population before it is lost to the star by gas drag. This concern is ameliorated by enhancements in particle density (pileups) that may occur as particles drift radially inward (Youdin & Shu 2002; Youdin & Chiang 2004), and by the reduction of drift speeds brought about by multiple particle sizes (Bai & Stone 2010a, see their Figure 8).

Regardless of which scenario nature prefers—particle concentration by the SI; dissipative gravitation into rings; or dynamical collapse of a vertically settled sublayer, which is the subject of this paper—all the proposed ways of forming planetesimals depend on knowing how far down dust settles and what maximum dust-to-gas ratios μ₀ can be attained at the midplane. Our order-of-magnitude estimate in Equation (11) requires testing. Among the most realistic simulations of particle settling are those by Johansen et al. (2009) and Bai & Stone (2010a), both of which concentrated on the SI. Johansen et al. (2009) reported that super-centimeter-sized particles settled into sublayers in which the midplane-averaged μ₀ ranged from 0.6 to 9.0 as the bulk metallicity ranged from 1 to 3× solar. Bai & Stone (2010a) found that the highest τ_s particles settled the most, driving turbulence that lofted smaller τ_s particles to greater heights. They argued that in their simulations, all of which were characterized by max τ_s ⩾ 0.1, particles were so strongly stirred by the SI that the KHI never manifested.

To complement these studies, we would like to understand the settled equilibrium states of disks composed entirely of small particles, well but not perfectly coupled to gas (0 < max τ_s ≪ 1), isolated from the complicating effects of the SI but not other instabilities like the KHI. Some previous attempts in this regard relied on assumed forms for the density profile of settled dust. Barranco (2009) presumed the dust density profile was Gaussian in shape, and did not seek to determine the maximum value of μ₀ per se. Chiang (2008) and Paper I, following Sekiya (1998) and Youdin & Shu (2002), assumed the dust density profile had a spatially constant Richardson number. We found in Paper I that under this assumption, in a disk of bulk solar metallicity (Σ_d/Σ_g = 0.015), the sublayer could remain KH stable for μ₀ as high as 8—a value that is nearly an order of magnitude higher than our crude estimate in Equation (11), and as such lowers the hurdle to forming planetesimals by gravitational instability.

Weidenschilling (2006, 2010) did not assume a shape for the marginally stable density profile to which dust relaxes. Instead the profile was calculated from a one-dimensional (1D) model that balanced the downward flux of particles by gravity (including vertical self-gravity) with the upward flux due to turbulent diffusion. The model used prescriptions for KH turbulence, with a number of parameters chosen to match the 3D calculations of Cuzzi et al. (1993)—which also relied on a prescribed form of the turbulence. Weidenschilling (2006) found that the density of millimeter-sized particles could exceed the Toomre density at r = 3 AU in a disk approximately 1.3× as massive as our minimum-mass solar nebula, and whose bulk metallicity Σ_d/Σ_g = 0.054 ≈ 3.6× solar. See his Figure 10, but note that 3.6× solar metallicity ≈16× his "nominal" abundance of solids, and that his "critical" density is 3× larger than our Toomre density.

Like Weidenschilling (2006, 2010), we seek the marginally stable state to which dust settles, in the limit that particles are well but not perfectly coupled to gas (0 < τ_s ≪ 1). We improve upon these earlier studies by not prescribing or parameterizing the turbulence, but by letting turbulence arise and evolve naturally from our 3D integrations of the standard fluid equations. In particular, our work is free of the popular but untested assumption that dust settles until the Richardson number equals a constant everywhere. We allow dust grains to fall until they are stopped by whatever instabilities they self-generate. In the calculations presented here we assume that dust begins well mixed with gas in a Gaussian density profile, and then follow the dust into whatever non-Gaussian distribution it seeks to relax. Although we try only a Gaussian initial profile, our method accommodates arbitrary initial conditions.

At the heart of our approach lie two codes. The first code is in one dimension (z) and computes the vertical drift of dust grains at their terminal velocities. Though incapable of deciding whether the density profiles it generates are prone to the KHI (or any other instability), the 1D code can evolve dust profiles for the entire settling time t_settle, which can be arbitrarily long for arbitrarily small grains. The task of assessing stability is reserved for the second code: the spectral, anelastic, shearing box code of Barranco (2009) which treats gas and dust in the perfectly coupled τ_s = 0 limit. Though incapable of allowing dust to settle out of gas, the 3D code accounts for the complicated interplay of vertical shearing and rotational effects to decide whether a given dust layer overturns from the KHI (or some other instability). It tests dynamical stability by running for dozens of dynamical times t_dyn = Ω⁻¹_K. Our procedure involves alternating between these two codes: allowing dust to settle over some fraction of the settling timescale t_settle using the 1D code; passing the results of the 1D code to the 3D code and allowing the dust profile to relax dynamically over timescales t_dyn; passing the results of the 3D code back to the 1D code for further sedimentation on the settling timescale; and so on, back and forth, until the midplane dust-to-gas ratio stops increasing, at which point the marginally stable state is identified.

In Section 2 we describe our method in full. Results are presented in Section 3, extended in Section 4, and summarized and discussed in Section 5.

As sketched in Section 1, to find the marginally stable state to which small dust grains relax, we alternate between two codes: a 1D code that tracks how dust drifts toward the midplane on the settling timescale t_settle and a 3D shearing box code developed by Barranco (2009) that allows dusty gas to stabilize on the dynamical timescale t_dyn ≪ t_settle. The 3D code integrates the anelastic fluid equations for perfectly coupled dust and gas using a spectral method. It includes a background radial pressure gradient to drive a vertical shear. Details about the 3D code are in Barranco & Marcus (2006), Barranco (2009), and Paper I.

Dust and gas are initially well mixed with a spatially constant density ratio: [ρ_d(z)/ρ_g(z)]_init ≡ μ_init = constant. We set μ_init equal to either solar metallicity (μ_init = 0.015; Lodders 2003; see also footnote 7) or four times solar metallicity (μ_init = 0.06). To determine the initial form of the dust density profile ρ_d(z), we solve the equation for vertical hydrostatic equilibrium where gas is assumed to be initially isothermal:

$$\begin{equation} \frac{c_{\rm s}^2}{\rho _{\rm g}+\rho _{\rm d}}\frac{\partial \rho _{\rm g}}{\partial z} = -\Omega ^2_{\rm K}z, \end{equation} \tag{ 12 }$$

whence

$$\begin{equation} \rho _{\rm d}= \mu _{\rm init}\rho _{\rm g0} \exp \left[-\frac{(1+\mu _{\rm init})z^2}{2H^2_{\rm g}}\right] \end{equation} \tag{ 13 }$$

for constants μ_init, a characteristic initial height H_g ≡ c_s/Ω_K, and the midplane gas density ρ_g0.

Taking μ_init to be constant over all gas scale heights is a simplifying but probably unrealistic assumption. Even if particles of fixed size s carry the bulk of the disk's solid mass, such particles would not likely begin well mixed with gas everywhere. Densities at altitude may be too low to permit particles of size s to coagulate, and to keep such particles aloft once formed. Nevertheless, the error accrued by our assumption of constant μ_init is small insofar as increasingly small amounts of mass are contained at larger altitudes. Moreover, we will find evidence that the final marginally stable state to which dust relaxes is insensitive to initial conditions (Section 5.1). In any case, we will point out in Section 3 which features of the evolving dust profile are artifacts of our assumed initial condition.

Equation (13) defines the initial dust profile used by the 1D code, whose grid extends from z = 0 to z = 3H_g. The 1D Lagrangian code uses particles to track the motion of dust mass. Each particle represents the same amount of dust mass. Any dust density profile ρ_d(z) can be converted into particle positions and back again. The closer particles are spaced, the greater is ρ_d.

Starting with Equation (13), we proceed as follows.

• 1.  
  1D code: Initialize positions of dust particles and establish hydrostatic equilibrium for the gas.Given ρ_d(z), calculate the positions of ∼60,000 particles in the 1D code. Also determine the hydrostatic gas density ρ_g(z) at each particle's position by solving Equation (12).⁹
• 2.  
  1D code: Settle dust particles by one timestep Δt. By equating the vertical gravitational force ∝z to the Epstein drag force F_D ∝ ρ_gv_rel (e.g., Weidenschilling 1977), assign terminal velocities

  $$\begin{equation} v_{\rm rel} \propto \frac{z}{\rho _{\rm g}} \end{equation} \tag{ 14 }$$

  to each particle. Advect each particle vertically downward by a distance v_relΔt, where Δt is chosen small enough that particles do not overtake one another. Note that the coefficient of proportionality on the right-hand side of Equation (14), which depends on quantities such as the disk mass parameter F and grain properties s and ρ_s, does not affect the shapes of the density profiles generated so long as it is the same for all particles.Bin particle positions and recalculate ρ_d(z).
• 3.  
  1D code: Repeat steps (1) + (2) until the midplane dust density ρ_d(0) rises by 30%.
• 4.  
  1D → 3D code: Insert results of the 1D code into the 3D code and run the 3D code for 10 orbits. Let z_max be the position of the highest particle in the 1D code, and set the dimensions of the shearing box in the 3D code to be (L_r, L_ϕ, L_z) = (1.455, 2.91, 4)z_max, resolved by (N_r, N_ϕ, N_z) = (32, 64, 128) gridpoints.¹⁰ Initialize the 3D code by assigning μ(z), as calculated by the 1D code, to each horizontal gridpoint (r, ϕ).¹¹ The 3D code initializes the remaining variables—velocity, gas density, temperature, and enthalpy—to ensure dynamical equilibrium; see Section 2.2 of Paper I.The background radial pressure gradient is parameterized by the variable v_max:

  $$\begin{equation} -\frac{1}{\rho _{\rm g}}\frac{\partial P}{\partial r} = 2\Omega _{\rm K}v_{\rm max}\,. \end{equation} \tag{ 15 }$$

  We fix v_max = 0.025c_s for all simulations. Physically, v_max ∼ c²_s/(Ω_Kr) represents the difference in azimuthal velocity between pressure-supported dust-free gas and a strictly Keplerian flow.¹²Before running the 3D code, perturb μ(r, ϕ, z) by an amount

  $$\begin{eqnarray} \Delta \mu (r,\phi,z) &=& A(r,\phi) \mu (z) [ \cos (\pi z / 2 z_{\rm max})\nonumber\\ && + \sin (\pi z/2 z_{\rm max}) ], \end{eqnarray} \tag{ 16 }$$

  where A(r, ϕ) is a random variable constructed in Fourier space (see the discussion following Equation (31) of Paper I). Fix the root mean square of the perturbations to be A_rms ≡ 〈A²〉^1/2 = 10⁻³.Run the 3D code for 10 orbits.
• 5.  
  3D code: Assess stability. Extend simulations beyond 10 orbits as necessary to make this assessment. Label the dust profile "KH-unstable" if the horizontally averaged dust-to-gas ratio at the midplane

  $$\begin{equation*} \langle \mu (z=0) \rangle \ \rm{as \ a \ function \ of}\ {\it t} \end{equation*}$$

  decreases by more than 15%. Otherwise, monitor the horizontally averaged vertical kinetic energy at the midplane:

  $$\begin{equation*} \langle \mu v^2_z (z=0)\rangle /2\ \rm{as \ a \ function\ of}\ {\it t}. \end{equation*}$$

  If 〈μv²_z〉/2 monotonically decreases or levels off, label the dust profile "KH-stable." If 〈μv²_z〉/2 is increasing toward the end of the simulation, extend the integration an additional 10 orbits and re-assess stability. Repeat step (5) as necessary.
• 6.  
  If "KH-unstable," stop. Identify the last KH-stable dust profile, generated in the iteration just previous to that of the KH-unstable simulation, as the "marginally stable state."¹³
• 7.  
  If "KH-stable," pass results of the 3D code back to the 1D code and return to step (1). Fit a polynomial μ_poly(z) to the final, horizontally averaged dust-to-gas ratio 〈μ(z)〉 as calculated by the 3D code. Adjust the order of the polynomial to capture all features of the profile. If one polynomial is insufficient, use two to create a piecewise function. Convert μ_poly(z) to ρ_d(z) by assuming the gas profile to be Gaussian (see footnote 9): ρ_d(z) = μ_poly(z) · ρ_g0exp [ − z²/(2H²_g)]. Using this ρ_d(z), return to step (1) for the next iteration.

Our method and results apply to any location in a disk of any mass (arbitrary r and F), provided our input assumptions that self-gravity is negligible and v_max/c_s = 0.025 are satisfied. They also apply to any particle size to the extent that the disk's solid mass is concentrated in particles of a single size (so that z and ρ_g uniquely determine v_rel; Equation (14)), and to the extent that such particles are undisturbed by SIs (Sections 1 and 5.3). Another way of saying all this is to note that our calculations are carried out in dimensionless units.

To orient the reader, in Figure 1 we show results obtained from the 1D code only. The dust is initially well mixed with gas at solar metallicity (μ_init = 0.015). As dust settles and the midplane dust-to-gas ratio μ₀ increases, sharp cusps appear at the edges of the dust layer where particles pile up vertically. Pileups occur because particle fluxes ρ_d|v_rel| ∝ μ|z| increase with increasing height |z|. This follows from our assumption that μ_init is constant, which as noted at the beginning of Section 2 may not be realistic.

Zoom In Zoom Out Reset image size

Unlike us, Garaud & Lin (2004) did not find vertical pileups at the edges of their layer because they chose their initial dust profile to have a scale height equal to 0.1H_g. Their initial μ profile decreased with |z| more quickly than 1/|z|, and thus did not satisfy the condition for pileups. We verified this by inserting their initial profile into our 1D code.

The shapes of the settled dust profiles μ(z) and their relative spacing in time are independent of the dust internal density ρ_s, dust particle size s, and the scaling parameter F for disk mass. Changing these parameters only alters the absolute physical time elapsed (Equation (3)). Relative time is tracked by the dimensionless parameter f ≡ t/t_settle, labeled on this and many subsequent figures.

Below we compare these 1D-only results to those that include the full 3D dynamics. The solar metallicity case is described in Section 3.1. The metal-rich case (μ_init = 0.06) is presented in Section 3.2.

3.1. Solar Metallicity

Figure 2 traces the evolution of dust that starts well mixed with gas at solar metallicity. Plotted are several KH-stable curves from the 3D code resulting from step (5) of our procedure. For ease of comparison with the purely 1D results, the relative timestamps in Figure 2, measured by f, coincide with those in Figure 1. The leftmost curve at f = 1.0 represents the marginally stable state identified using our standard procedure. This state achieves a midplane dust-to-gas ratio of μ₀ = 2.45, about an order of magnitude below the value required for gravitational instability (Equation (10)). In Section 4, we extend our procedure to see if we might achieve still higher dust-to-gas ratios.

Zoom In Zoom Out Reset image size

Comparing Figures 1 and 2, we see that the (possibly unrealistic) pileups at the edges of the dust layer do not survive in the dynamical 3D code. By f ≈ 0.44, the pileups are nearly gone. At this point, the vertical extent of the dust layer z_max has shrunk to ∼0.1H_g, and the only pileup present is the one at the midplane.

The instability that eliminates the pileups at the edges of the layer is likely related to the Rayleigh–Taylor instability (RTI), triggered by heavy fluid lying on top of lighter fluid, and we will refer to it henceforth as such. The RTI originates locally at the edges of the dust layer. By contrast, the midplane is relatively stable (at least until the marginally stable state is reached). Another way of seeing this is to note that midplane dust-to-gas ratios μ₀ in Figures 1 and 2 agree to within 25%. Closer examination reveals that those in Figure 2 are consistently higher. This suggests that the RTI transfers some of the dust in the pileups to the midplane.

Figure 3 confirms this transfer mechanism. The top middle panel shows that over the course of a 20-orbit-long 3D simulation (iteration #6, occurring at a time f = 0.31, out of a total of 19 iterations), dust is redistributed from the layer's edges to the midplane, raising μ₀ by about 20%. Note that the effect of the RTI has been to transport dust toward the midplane, not to higher altitudes. The RTI is confined to where dust is unstably stratified (increasing total density in the direction opposite to gravity).

Zoom In Zoom Out Reset image size

Compare this behavior with that in the top row of Figure 4, which documents a later iteration, #16. The top middle panel shows that an instability has occurred near the edges of the dust layer. Dust is redistributed to higher, not lower, altitudes. The midplane is not affected. The instability at this relatively late stage of settling is probably driven by the vertical shear associated with strong density gradients at the edges of the layer, and we will refer to it henceforth as the KHI. As a result of the KHI, gradients in density and velocity are reduced.

Zoom In Zoom Out Reset image size

The marginally stable state identified using our standard procedure is displayed in Figure 5. The bottom panels show that during the last iteration #19, the usual 30% increase in the midplane μ₀ (left bottom) results in a KH-unstable profile (middle bottom). In the top panels, we re-do iteration #19, this time incrementing μ₀ by only 10% (left top). The resultant profile is KH stable (middle and right top panels), and has 〈μ₀〉 = 2.45. In Section 4.2, we modify our standard procedure and extend it to later times to achieve still higher dust-to-gas ratios in stable flows.

Zoom In Zoom Out Reset image size

The μ-profiles in Figures 2–5 betray oscillations just inside the edges of the dust layer. We believe these ripples are artificial because when each first appears, it spans only a few grid points of the 3D code: see the f = 0.054 profile of Figure 2, which shows two nascent ripples. The features probably arise because the truncated Chebyshev series used to model the flow in z has too few terms to adequately capture the steep vertical density gradient (Gibbs 1898). Originating in the 3D code, the ripples are then amplified as mini-pileups in the 1D code. We could have tried to smooth away these oscillations by reducing the order of our polynomial fit (step 7 of our procedure), but chose instead to retain all features of the dust profile generated by both codes to minimize bias. In any case the oscillations are eventually erased by instabilities during the later stages of settling (Figure 2). In and of themselves the oscillations do not appear to introduce instabilities, which as discussed above are triggered instead by smooth density gradients—realistically computed—at the boundaries of the layer (top rows of Figures 3 and 4).

3.2. Metal-rich Case: 4 × Solar Metallicity

Figure 6 follows the evolution of dust that is initially well mixed with gas at 4× solar metallicity. It shares the same timeline as Figures 1 and 2. Thus, the last profile marked f = 1.1 in Figure 6 is attained at a time 10% later than that marked f = 1.0 in the other figures. This last profile is the marginally stable state identified using our standard procedure, for the case of supersolar metallicity. It achieves a midplane dust-to-gas ratio of μ₀ = 20.3—large enough to exceed the Toomre threshold in a disk that has twice the gas content of the minimum-mass solar nebula (F = 2 in Equation (10)).

Zoom In Zoom Out Reset image size

The evolution of the metal-rich disk over the course of a total of 21 iterations—some of which are sampled in Figures 7–9—is similar to that of the solar metallicity disk, with two notable differences. When the unstably stratified pileups of dust collapse (iteration #4, shown in the top row of Figure 7), enough dust is transferred to the midplane that μ attains an appreciable maximum there. This bump contrasts with the nearly flat profile seen for the solar metallicity run (Figure 3), and persists at least through iteration #16 (Figure 8). A second difference is that in every KH-stable simulation following iteration #13, the vertical kinetic energy, although it eventually levels off, ends orders of magnitude higher than where it began (Figure 8, and the top row of Figure 9). Some currents and/or turbulence appear to be sustained as the state of marginal stability is approached.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Related to this second point, we should acknowledge that our standard procedure ignores whatever velocities are present at the end of a given 3D simulation when initializing the velocities of the subsequent 3D simulation. That is, with every iteration, velocities are set anew according to Equation (5), with vertical and radial velocities reset to zero. The assumption we make in our standard procedure is that whatever velocities are maintained in a KH-stable layer do not stop dust from settling at the local terminal velocity v_rel. A crude attempt at relaxing this assumption is made in Section 4.2.

In Section 4.1 we test the robustness of our results against the size of our simulation box. In Section 4.2 we modify the procedure of Section 2, pushing to still higher dust-to-gas ratios at the midplane and revising our identification of marginally stable states.

4.1. Bigger Box Runs

4.2. Refining the Marginally Stable State Using a Modified Settling Procedure

Using our standard procedure of Section 2, we can only provisionally identify marginally stable dust profiles. The identification is provisional because by settling all dust particles at their local terminal velocities v_rel, we wind up with edges so unstable that they also destabilize the midplane. In reality, dust particles at the edge may stop settling because they attain a state of marginal stability first, remaining lofted up by the gas motions they stir up locally, and leaving dust particles near the midplane free to keep settling. In other words, marginal stability may be reached sequentially, starting from the edges and ending at the midplane. Our standard procedure does not allow for this kind of gradual evolution because the 1D code settles all dust particles at their local terminal speeds regardless of their location. In this sense our standard procedure is too blunt because it does not allow for slower settling at the stirred-up edges and faster settling at the quiescent midplane. True marginally stable states should have midplane dust-to-gas ratios even higher than the maximum ones displayed in Figures 2, 5, 6, and 9.¹⁴

To remedy this shortcoming, we modify our procedure by applying a weighting function 0 < W(z) ⩽ 1 to each dust particle's settling velocity. Starting with a KH-stable state near the end of our standard sequence of iterations, we continue to let dust particles drift to the midplane in the 1D code, not at v_rel(z) but at the weighted velocity W(z)v_rel(z). We use either a Fermi function

$$\begin{equation} W (z) = \frac{1}{1 + \exp [(z-z_{\rm 50})/z_{\rm w}]}, \end{equation} \tag{ 17 }$$

described by two parameters z₅₀ and z_w, or a Gaussian

$$\begin{equation} W (z) = \exp \big(-z^2 /2z^2_{\rm w}\big), \end{equation} \tag{ 18 }$$

described by a single parameter z_w. The choice of weighting function is somewhat arbitrary; it depends on the shape of the dust profile to be settled and is made case-by-case according to considerations outlined below. The intent of the weighting function is to slow the steepening of density gradients at the dust layer's edges and thereby prevent the edges from destabilizing the entire layer.

We start with the KH-stable profile in iteration #18 of our solar metallicity sequence (black solid curve in the top left panel of Figure 10). The dust layer is characterized by a "core" from z = 0 to z ≈ 0.5z_max over which 〈μ〉 is fairly constant, and an "edge" from z ≈ 0.5z_max to z = z_max over which the dust content drops to zero. Because the instabilities that threaten to disrupt the layer originate in the edge and not the core, we seek a weighting function W(z) that slows the downward drift of dust in the edge but not in the core. At the same time W(z) should not be so strongly weighted toward the midplane that the core disconnects from the edge and opens a local gap in dust content. We find upon experimentation that a Gaussian does not have enough flexibility to meet these requirements for this particular iteration. However, a Fermi function with z₅₀ = 0.005H_g—corresponding approximately to the boundary of the core—and z_w = 0.05z₅₀ seems to work well (blue dashed curve in the top left panel of Figure 10). We use this weighting function to settle the dust profile until its midplane μ₀ increases by 30% to a value of 2.9 (red dot-dashed curve). This settled layer remains KH stable (top middle and right panels of Figure 10)—unlike the layer settled without the weighting function (bottom panels of Figure 5).

Zoom In Zoom Out Reset image size

In the new profile settled with our modified procedure, the distinction between the core and the edge is no longer so sharp. Thus, to settle this new profile even further, a simple Gaussian for the weighting function seems to suffice. Choosing z_w = 0.0025H_g ≈ 0.25z_max, we attempt to increase the midplane μ₀ yet again by 30%, but find the resultant profile to be KH unstable (bottom panels of Figure 10).

Similar results are obtained for the metal-rich case (Figure 11). Using Gaussian weighting functions we are able to push the midplane dust-to-gas ratio μ₀ to a new record of 26.4, which is 30% greater than the value attained using our unweighted standard procedure.

Zoom In Zoom Out Reset image size

To form rocky planets and gas giant cores, dust must amass in a circumstellar disk. In the classic scenario for forming planetesimals, dust settles vertically toward the midplane into an ever thinner and denser layer that eventually becomes gravitationally unstable (Safronov 1969; Goldreich & Ward 1973). Toomre's criterion for gravitational instability (GI) is satisfied for midplane dust-to-gas ratios (ρ_d/ρ_g)₀ ≡ μ₀ ≳ μ_0,Toomre, where μ_0,Toomre ≈ 34 for a minimum-mass nebula orbiting a solar-mass star (equation 10; note that μ_0,Toomre is nearly independent of disk radius). For comparison, in a disk of well-mixed dust and gas at solar abundance, μ₀ ≈ 0.015 (Lodders 2003). Whether dust can accumulate until its density increases by more than three orders of magnitude depends on how turbulent the ambient gas is. Even supposing that gas in certain regions of the disk is not intrinsically turbulent (e.g., because it is too weakly ionized to support the magnetorotational instability), the dust itself may excite turbulence in gas by a Kelvin-Helholtz-type/shearing instability (KHI). The KHI is triggered when the velocity gradient between dust-rich gas at the midplane and dust-poor gas at altitude becomes too large. Barring gravitational instability, dust should settle to a state that is marginally stable against the KHI. The question of whether gravitational instability is viable translates into the question of whether the state that is marginally stable to the KHI has μ₀ ≳ μ_0,Toomre (this is a necessary but not sufficient criterion for the formation of planetesimals by gravitational collapse; see footnote 5).

In this paper, we sought out such marginally stable states by numerical simulation. Starting with dust well mixed with gas at either bulk solar or supersolar metallicity, we allowed dust to settle vertically until dynamical instabilities prevented the midplane density from increasing further. We tracked the approach to the marginally stable state using a combination of a 1D settling code and a 3D shearing box code, working in the limit that particles are small enough not to excite streaming instabilities (SIs). All the instabilities that afflicted our dust layer originated at the layer's edges, where dust density gradients were steepest. We found evidence for two kinds of instabilities: the usual KHI driven by vertical shear, and the Rayleigh-Taylor instability driven by the weight of piled-up dust. These instabilities were mostly confined to the dust layer's top and bottom surfaces, leaving dust near the midplane free to settle but occasionally speeding up the accumulation of solids by transferring dust from pileups downward. The midplane density stopped increasing when the dust layer became so thin that instabilities at the edges threatened to overturn the entire layer.

Using our standard procedure of Section 2, we attained maximum dust-to-gas ratios of μ₀ ≈ 2.45 and 20.3 for the cases of solar and 4× solar bulk metallicity, respectively (Figures 5 and 9). These values are lower limits because in our standard procedure dust particles at the layer's top and bottom faces keep settling until they excite instabilities so vigorous that dust at the midplane is stirred up. In reality, dust at the layer's edges may reach a state of marginal stability and stop settling, leaving dust near the midplane free to settle further. We modified our procedure in Section 4.2 to try to account for this effect, reaching μ₀ ≈ 2.9 and 26.4 for the two metallicity cases (Figures 10 and 11). These values are still lower limits because our simulations omit self-gravity. But the correction for self-gravity should be small for the solar metallicity disk, on the order of 10% (∼2.9/34). For our supersolar metallicity disk, the correction for self-gravity might be on the order of unity (∼26.4/34)—although it might also be much higher, as Sekiya (1998) and Youdin & Shu (2002) showed that vertical self-gravity can yield a singularity in μ₀.

We conclude that a minimum-mass disk of bulk (height-integrated) solar metallicity orbiting a solar-mass star cannot form planetesimals by self-gravity alone: even neglecting turbulence intrinsic to gas, the KHI would force the midplane dust density to fall short of the Toomre density by about an order of magnitude. Our results make clear what changes to the circumstellar environment would be needed for self-gravity to prevail. To attain the Toomre density in a minimum-mass gas disk, the bulk metallicity would need to be enhanced over solar by a factor of a few ≲4. For disks with total mass (gas plus dust) greater than that of the minimum-mass solar nebula, the required degree of metal enrichment would be lower.

Our results agree with those of the prescriptive model of Weidenschilling (2006), who found that the density of millimeter-sized particles (τ_s ∼ 0.001) at r = 3 AU in a disk for which ρ_g = 1.6 × 10⁻⁹ g cm⁻³ (F ≈ 1.3) and Σ_d/Σ_g ≈ 0.015 (solar metallicity) fell short of the Toomre density by about a factor of 10. When the bulk metallicity Σ_d/Σ_g increased to 0.054, the Toomre density was exceeded by a factor of 3.

5.1. How Spatially Constant is the Richardson Number?

In Paper I, as in previous works (Sekiya 1998; Youdin & Shu 2002; Youdin & Chiang 2004), all dust profiles were assumed to have spatially constant Richardson numbers Ri. The dust profiles we have computed are free of this assumption, whose validity we can now test.

We calculate Ri(z) for our marginally stable states, derived under both standard (Section 2) and modified (Section 4.2) procedures. To compute the numerator (Brunt frequency) of Ri in Equation (4), we use the horizontally averaged dust-to-gas ratio 〈μ(z)〉, computing derivatives using centered differences and assuming the gas to obey a Gaussian density profile (see footnote 9). To compute the denominator (vertical shearing rate) of Ri, we also use 〈μ(z)〉, inserting it into Equation (5) and computing therefrom the velocity derivative. Of course we could also compute the denominator more directly by using the simulation output itself for v_ϕ—this alternative approach turns out to give identical results for the solar metallicity disk, but for the metal-rich disk the Ri(z) so generated varies much more erratically. As noted in Section 3.2, the metal-rich disk sustains gas motions well above those we put in as seed noise. These motions are not strong enough to overturn the dust layer but they are large enough to render Ri highly variable, both in time and space. By not using the simulation data for v_ϕ(z) and relying instead on the better behaved 〈μ(z)〉, we effectively average Ri in time and space.

Results for the solar metallicity runs are shown in Figure 12. We plot Ri only where 〈μ〉-gradients are large enough to be computed reliably—thus we avoid regions closest to the midplane. Although we find that Ri is not a strict constant, it varies only between 0.1 and 0.3 within a large fraction of the edges of the dust layer—precisely where instabilities, presumably shear-driven, have rearranged dust. This result compares favorably with Paper I, where we found that a solar metallicity disk has a critical Richardson number of Ri_crit ≈ 0.2.

Zoom In Zoom Out Reset image size

Evidence for a constant Ri within the edges of the dust layer is even stronger for the metal-rich disk, as shown in Figure 13. Here, Ri hovers near 0.5 over much of the edges—again consistent with Paper I. See Figure 6 of that paper; admittedly the curve for Ri_crit(Σ_d/Σ_g) in Paper I needs to be extrapolated to the supersolar metallicity considered here, Σ_d/Σ_g = 0.06.

Zoom In Zoom Out Reset image size

The Ri(z) profiles in Figures 12 and 13 differ from those in Figure 3 of Bai & Stone (2010a); see the dashed curves corresponding to their 3D simulations, all of which include marginally aerodynamically coupled particles. These differences support their arguments that their simulations were not afflicted by the KHI.

In summary, the assumption made in other studies that well-coupled particles settle into a layer for which Ri is spatially constant does not appear too bad. This is welcome news, not least because it implies that the final marginally stable states are robust against uncertainties in initial conditions (Section 2). One caveat is that we have not tested this assumption for those regions closest to the midplane, as they could not relax by our method before being disrupted by instabilities at the edges. This is an area where more work can be done; see Section 5.3. Another caveat, supported independently by Paper I, is that the critical value of Ri to which dust relaxes is not unique, but increases with bulk metallicity Σ_d/Σ_g. For a solar metallicity disk, Ri_crit ≈ 0.2, but for one having 4× solar metallicity, Ri_crit ≈ 0.5. This trend has not yet been explained.

5.2. The Super-linear Relation Between Maximum Dust-to-gas Ratio μ₀ and Bulk Metallicity Σ_d/Σ_g

The degree to which Ri is constant is related to the scaling between the maximum midplane dust-to-gas ratio μ₀ and bulk metallicity Σ_d/Σ_g. Naively it might be expected that μ₀ scales linearly with Σ_d/Σ_g—a greater total amount of metals simply yields a proportionately dustier midplane—and indeed a linear relation is implied by our order-of-magnitude estimate in Equation (11). But we did not find a linear trend in our simulations. We find instead that the relation is super-linear: a factor of 4 increase in Σ_d/Σ_g results in a factor of 9.1 increase in maximum μ₀ (26.4 versus 2.9).

A super-linear trend is predicted by theories assuming a constant Richardson number (Sekiya 1998; Youdin & Shu 2002). The large gain in midplane density afforded by a comparatively modest increase in bulk metallicity is partly what makes planetesimal formation by gravitational instability so attractive. Increasing Σ_d/Σ_g does more than just increase the total amount of metals in the disk—it also helps to stabilize it, by decreasing the vertical shear ∂v_ϕ/∂z. In the limit μ₀ ∼ (Σ_d/Σ_g)H_g/Δz ≫ 1, where Δz is the dust layer thickness, we have from Equation (5)

$$$ \frac{\partial v_\phi }{\partial z} \sim \frac{\eta \Omega _{\rm K}r / \mu _0}{\Delta z} \sim \frac{\eta \Omega _{\rm K}r}{H_{\rm g}} \frac{1}{\Sigma _{\rm d}/\Sigma _{\rm g}} $$$

which decreases with increasing Σ_d/Σ_g. By comparison, the Brunt frequency [(g/ρ)∂ρ/∂z]^1/2 ∼ [(Ω²_KΔz/μ)μ/Δz]^1/2 ∼ Ω_K hardly changes with Σ_d/Σ_g. Thus, the Richardson number increases as Σ_d/Σ_g increases; the enhanced stability allows Δz to decrease; and thus μ₀ ∝ Σ_d/Δz scales super-linearly with Σ_d.

The above order-of-magnitude relations show qualitatively how a super-linear trend follows from the decrease in vertical shear, and the consequent increase in stability, brought about by an increase in bulk metallicity. However, these simple relations are not enough to quantify the super-linear trend because Δz appears to have cancelled out of both ∂v_ϕ/∂z and the Brunt frequency. This difficulty is avoided in a more formal derivation of the relation between μ₀ and Σ_d/Σ_g, made under the assumption of constant Ri, as described in Appendix B.

We note further that μ₀ scales with the inverse of the radial pressure gradient parameter v_max/c_s (equivalently η) in the same super-linear way as for Σ_d/Σ_g. The smaller is v_max/c_s, the greater is the maximum μ₀ attainable; the relation between these quantities is derived under the assumption of constant Ri in Appendix B. Thus, we expect our numerical results for max μ₀ (2.9, 26.4) to depend sensitively on our assumed value of v_max/c_s = 0.025. (Bai & Stone 2010c also reported that the degree of clumping caused by the SI increased strongly with decreasing v_max/c_s.)

5.3. Future Directions

5.4. Connection to Observations and the Need for Supersolar Metallicity

We thank Xue-Ning Bai, John Johnson, Eve Ostriker, Jim Stone, and Neal Turner for discussions, and Tristan Guillot for the data in Table 1. Xue-Ning Bai, Anders Johansen, Jim Stone, and Andrew Youdin provided valuable feedback on a draft version of this paper. We are grateful to Stuart Weidenschilling for an insightful referee's report that put our work into better context. This research was supported by the National Science Foundation, in part through TeraGrid resources provided by Purdue University under grant number TG-AST090079. A.T.L. acknowledges support from an NSF Graduate Fellowship.

For numerical estimates in this paper, we adopt the standard disk model derived in the review by Chiang & Youdin (2010). The disk has surface densities

$$\begin{eqnarray} &&\Sigma _{\rm g} = 2200 \,F \left(\frac{r}{\rm AU} \right)^{-3/2} {\rm \,g}{\rm \,cm}^{-2} \end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray} &&\Sigma _{\rm d} = 33 \,F\, Z_{\rm rel}\left(\frac{r}{\rm AU} \right)^{-3/2} {\rm \,g}{\rm \,cm}^{-2} \end{eqnarray} \tag{ A2 }$$

in gas (g) and dust (d). The dimensionless parameters F and Z_rel ≡ (Σ/Σ_g)/0.015, typically of order unity, describe how much total mass the disk has relative to the minimum-mass solar nebula, and how metal-rich the disk is compared with a gas of solar abundances, respectively. The minimum-mass solar nebula (F = 1, Z_rel = 1) uses a condensate mass fraction for solar abundances of Σ_d/Σ_g = 0.015 (Lodders 2003). Values of Z_rel>1 correspond to supersolar metallicities Σ_d/Σ_g>0.015. Integrated to r = 100 AU, Equation (A1) yields a total disk mass of 0.03F M_☉.

At the disk midplane, the gas temperature, scale height, and density are given by

$$\begin{eqnarray} T &=& 120 \left(\frac{r}{\rm AU} \right)^{-3/7} {\rm \,K} \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} H_{\rm g} &=& 0.022 r \left(\frac{r}{\rm AU} \right)^{2/7} \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} \rho _{\rm g0} &=& 2.7 \times 10^{-9} F \left(\frac{r}{\rm AU} \right)^{-39/14} {\rm \,g}{\rm \,cm}^{-3} \,. \end{eqnarray} \tag{ A5 }$$

These are adapted from Chiang & Goldreich (1997), adjusted for a disk obeying (A1)–(A2), orbiting a pre-main-sequence star of mass M_* = 1 M_☉, radius R_* = 1.7 R_☉, and temperature T_* = 4350 K.

We derive μ₀ as a function of Σ_d/Σ_g under the assumption of a constant Ri. Some evidence supporting a constant Ri was found in our simulations (Section 5.1). The density profile for constant Ri is used in a number of papers (Sekiya 1998; Youdin & Shu 2002; Paper I) and we begin by repeating the result, neglecting self-gravity as we have throughout our paper. The dust-to-gas ratio is given by

$$\begin{equation} \mu (z) = \left[\frac{1}{1/(1+\mu _0)^2 + (z/z_{\rm d})^2}\right]^{1/2} - 1, \end{equation} \tag{ B1 }$$

where

$$\begin{equation} z_{\rm d} \equiv \frac{Ri^{1/2}\,v_{\rm max}}{\Omega _{\rm K}} \end{equation} \tag{ B2 }$$

is a characteristic dust height and v_max = ηΩ_Kr (see Equations (5) and (6)) is a constant equal to the difference in azimuthal velocity between a strictly Keplerian flow and dust-free gas. The dust density drops to zero at

$$\begin{equation} z = \pm z_{\rm max} = \pm \frac{\sqrt{\mu _0(2+\mu _0)}}{1+\mu _0} z_{\rm d} \,. \end{equation} \tag{ B3 }$$

A comment on Equation (B1), in the limit that μ₀ ≫ 1: except where μ is nearly constant near z ≪ z_max/μ₀ and where it falls to zero near z = z_max, the shape of μ(z) is that of 1/z. This form follows simply from the constancy of Ri. Because the numerator of Ri is approximately constant with μ (Section 5.2), the denominator must be as well: ∂v_ϕ/∂z ∼ (v_max/μ)/z∼ constant, which implies μ ∝ 1/z. From this we can deduce the super-linear trend between μ₀ and Σ_d/Σ_g as follows. The integral of μ with respect to z is proportional to the total surface density of dust Σ_d. Because μ ∝ 1/z, flattening off as z decreases below z_max/μ₀, this integral varies as log μ₀. Then μ₀ ∝ exp Σ_d, crudely.

More formally, we have

$$\begin{equation} \Sigma _{\rm d}= 2 \int _0^{z_{\rm max}} \rho _{\rm d}\, dz = 2\rho _{\rm g0} \int _0^{z_{\rm max}} \mu \, dz \end{equation} \tag{ B4 }$$

where ρ_g0 is the midplane gas density, assumed constant because z_max ≪ H_g. The gas density profile always well approximates a Gaussian (see footnote 9), from which it follows that $\Sigma _{\rm g}\approx \sqrt{2\pi } \rho _{\rm g0} H_{\rm g}$. Then

$$\begin{equation} \frac{\Sigma _{\rm d}}{\Sigma _{\rm g}} = \sqrt{ \frac{2}{\pi } } \frac{1}{H_{\rm g}} \int _0^{z_{\rm max}} \mu \, dz \,. \end{equation} \tag{ B5 }$$

Inserting Equation (B1) into (B5), we have

$$\begin{equation} \sqrt{ \frac{\pi }{2} } \frac{H_{\rm g}}{z_{\rm d}} \frac{\Sigma _{\rm d}}{\Sigma _{\rm g}} = \log \big[1 + \mu _0 + \mu _0^{1/2}(2+\mu _0)^{1/2}\big] - \frac{\mu _0^{1/2}(2+\mu _0)^{1/2}}{(1+\mu _0)} \,. \end{equation} \tag{ B6 }$$

In the limit μ₀ ≫ 1, the exponential dependence of μ₀ on Σ_d/Σ_g is evident. Equation (B6) is plotted in Figure 14, with Ri = 0.25 and v_max/c_s = 0.025. Overlaid is the same equation but with varying Ri = Ri_crit ≈ 0.25(μ₀/9), the relation we found in Paper I (see Figure 5 of that paper). The two data points representing the maximum μ₀ achieved in this paper are also plotted. The data track the variable Ri_crit(μ₀) curve much better than the constant Ri curve.

Zoom In Zoom Out Reset image size

Finally note that H_g/z_d ∝ c_s/v_max enters into Equation (B6) the same way that Σ_d/Σ_g does. Thus, μ₀ increases super-linearly with c_s/v_max as well. This result leads us to suspect that our numerical results for μ₀ (2.9, 26.4) depend sensitively on our choice for v_max/c_s = 0.025. In this paper, we did not run simulations with different v_max/c_s and so did not test this suspicion.